Foo

Точкови оценки

Ако имаме наблюдения \(x_1, ..., x_n\) над случайна величина \(X \sim \mathcal{N}(\mu, \sigma^2)\), то можем да мислим на всяко наблюдение \(x_i\) като на наблюдение от случайна величина \(X_i\) разпределена като \(X\) и \(X_i\) са независими. Тоест, имаме независими и еднакво разпределени (н.е.р.) случайни величини \(X_1, ..., X_n \overset{distr}{\sim} X \sim \mathcal{N}(\mu, \sigma^2)\), като сме получили по едно наблюдение от всяка. Векторът \((X_1,...,X_n)\) наричаме случайна извадка.

Тогава \(\bar{x} := \frac{1}{n}\sum_{i=1}^n{x_i}\) е наблюдение от случайната величина \(\bar{X} := \frac{1}{n}\sum_{i=1}^n{X_i} \sim \mathcal{N}(\mu, \frac{\sigma^2}{n})\), понеже \[ E[\bar{X}] = E[\frac{1}{n}\sum_{i=1}^n{X_i}]=\frac{1}{n}\sum_{i=1}^n{E[X_i]}=\frac{1}{n}n\mu = \mu\\ Var[\bar{X}]=Var[\frac{1}{n}\sum_{i=1}^n{X_i}] = \frac{1}{n^2}n\sigma^2=\frac{\sigma^2}{n} \]

Забелязваме, че \(Var[\bar{X}] \xrightarrow[n\to \infty]{} 0\). Тоест, колкото повече наблюдения имаме, толкова повече \(\bar{x}\) има вероятност да е близо до истинския параметър \(\mu\), понеже ще е наблюдение от нормално разпределена случайна величина с дисперсия \(\frac{\sigma^2}{n}\) и очакване \(\mu\).

Можем да наблюдаваме този факт със следната симулация

# library("tidyverse")
# library("ggplot2")

# брой извадки които взимаме
n <- 1000

# 1000 пъти взимаме средно на 10, 100 и 1000 наблюдения от N(25, 40^2)
c(10, 100, 1000) %>%
  map(function(size) {
    # Взимаме средното на size на брой наблюдения от N(25, 40^2)
    map(1:n, ~mean(rnorm(size, mean = 25, sd = 40))) %>%
      unlist() %>%
      data.frame(x = ., size = as.character(size))
  }) %>%
  bind_rows() %>%
  ggplot() +
  geom_density(mapping = aes(x = x, fill = size)) +
  facet_wrap(~size) +
  ggtitle("Емпирично Разпределение на извадково средно при n=10, 100, 1000") +
  xlab("извадково средно") +
  ylab("емпирична плътност на извадково средно")

Още една визуализация

## Picking joint bandwidth of 0.712

\(\bar{X}\) наричаме точкова оценка за параметъра \(\mu\). Това, че \(Е(\bar{X}) = \mu\) наричаме неизместеност на оценката. Когато една точкова оценка \(f_{\theta}(\overrightarrow{X})\) приближава по вероятност своя параметър \(\theta\), тоест \[ \lim_{n\to \infty}Pr(|f_{\theta}(\overrightarrow{X}) - \theta| > \epsilon) = 0 \text{ за } \forall \epsilon > 0 \] казваме, че точковата оценка \(f_{\theta}(\overrightarrow{X})\) е асимптотично консистента.

Стандартното отклонение на оценката, обикновено наричаме стандартна грешка.

Интервални Оценки

Ако имаме наблюдения над случайна величина \(X \sim \mathcal{N}(\mu^?, \sigma^2)\) за която \(\mu^?\) е неизвестен параметър, но \(\sigma^2\) е известна (в практиката рядко ще знаем който и да е параметър), то знаем, че

\[ = := \frac{\bar{X} - \mu^?}{\sigma/\sqrt{n}} \sim \mathcal{N}(0, 1) \]

и следователно, можем да използваме определен квантил q от \(\mathcal{N}(0, 1)\), така че

\[ P(-q < \frac{\bar{X} - \mu^?}{\sigma/\sqrt{n}} < q) = \gamma \]

Така, с вероятност \(\gamma\) (ниво на достоверност) можем да сме сигурни, че

\[ \bar{X} - q \frac{\sigma}{\sqrt{n}} < \mu^? < \bar{X} + q \frac{\sigma}{\sqrt{n}} \]

В R за 95% интервал на достоверност, това ще стане по следния начин

n <- 1000
# сигма ни е известно, но мю не ни е и се опитваме да го оценим
known_sd <- 10
x <- rnorm(n, 14.55, known_sd)

emp_mean <- mean(x)
q <- qnorm(0.975)

c(emp_mean - q * known_sd / sqrt(n), emp_mean + q * known_sd / sqrt(n))
## [1] 13.88863 15.12822

В повечето случаи обаче, \(\sigma\) не ни е известно. Тогава, ако заместим с извадковото стандартно отклонение S, случайната величина \(\frac{\bar{X} - \mu}{S/\sqrt{n}}\) има \(\mathcal{T}\) разпределение с \(n - 1\) степени на свобода. В този случай, трябва просто да използваме квантилите от \(\mathcal{T}\) разпределение.

n <- 1000
x <- rnorm(n, 25, 7.654)

emp_mean <- mean(x)
emp_sd <- sd(x)

q <- qt(0.975, df = n - 1)

c(emp_mean - q * emp_sd / sqrt(n), emp_mean + q * emp_sd / sqrt(n))
## [1] 24.53077 25.47791

Или с t.test функцията

t.test(x)
## 
##  One Sample t-test
## 
## data:  x
## t = 103.61, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  24.53077 25.47791
## sample estimates:
## mean of x 
##  25.00434

Получената оценка за \(\mu\) се нарича интервална оценка. Ето как изглеждат интервалите на достоверност при Нормално и \(\mathcal{T}\) разпределение за случайната величина Z дефиниране по-горе.

# library("ggridges")

n <- 100000

list("normal" = rnorm(n), "t" = rt(n, 100000 - 1)) %>%
  imap(~data.frame(x = .x, distribution = .y)) %>%
  bind_rows() %>%
  ggplot(aes(x = x, y = distribution, fill = factor(stat(quantile)))) +
  stat_density_ridges(
    geom = "density_ridges_gradient",
    calc_ecdf = TRUE,
    quantiles = c(0.025, 0.975)
  ) +
  scale_fill_manual(
    name = "Probability", values = c("#FF0000A0", "#62b5b2", "#FF0000A0"),
    labels = c("(0, 0.025]", "(0.025, 0.975]", "(0.975, 1]")
  ) +
  xlab("Z") +
  theme_minimal()
## Picking joint bandwidth of 0.0898

Нека симулираме различни извадки от нормално разпределена случайна величина с \(\mu=3\) и \(\sigma^2 = 2\) и да видим в колко от случаите ще попадаме в интервала на достоверност.

simulate_t <- function(n, mu, sigma) {
  x <- rnorm(n, mu, sigma)
  b1 <- mean(x) - qt(0.975, df = n - 1) * (sd(x) / sqrt(n))
  b2 <- mean(x) + qt(0.975, df = n - 1) * (sd(x) / sqrt(n))
  between(mu, b1, b2)
}

# да преброим колко елементи са в интервала на достоверност
simulations <- replicate(10000, simulate_t(40, 4, 4))
mean(simulations)
## [1] 0.9532

Виждме, че наистина в около 95% от случаите попадаме в интервала на достоверност.

Тестване на хипотези върху една извадка

Хипотеза наричаме какво да е твърдение за функцията на разпределение \(F_X\) на случайната величина \(Х\) върху която имаме извадка \((X_1, ..., X_n)\).

Нека наблюдаваме данни идващи от нормално разпределение, примерно от \(\mathcal{N}(176, 6)\) за които не знаем истинската стойност на \(\mu\) (за момент забравяме, че сме ги генерирали с \(\mu=176\))

heights <- rnorm(100, 176, 6)

Ако примем, че истинското математическо очакване \(\mu\) е по-голямо или равно на 190, то случайната величина

\[ Z := \frac{\bar{X} - 190}{\sigma / \sqrt{n}} \] ще бъде нормално разпределена с \(\mu=0\) и \(\sigma^2=1\)

# z-statistic
z <- (mean(heights) - 190) / (6 / sqrt(100))
z
## [1] -20.99485

Можем да видим вероятността да имаме стойност z или по-малка използвайки функцията pnorm

pvalue <- pnorm(z)
# Вероятността да сме наблюдавали z или по-крайна стойност на z
pvalue
## [1] 3.654432e-98

Виждаме, че ако хипотезата, че \(\mu \ge 190\) e вярна, то е изключително малко вероятно да наблюдаваме \(z\)-статистиката която имаме (или по-крайна). Тогава, можем да отхвърлим хипотезата, че \(\mu \ge 190\). Но колко трябва да е малко нашето pvalue, така че да отхвърлим хипотезета \(H_0\) че \(\mu \ge 190\) и да приемем алтернативата \(H_1\) че \(\mu < 190\). Ако изберем някакво \(\alpha\), примерно \(\alpha = 0.05\) (ниво на значимост) можем да отхвърляме хипотезата \(H_0\) само когато \(pvalue < \alpha := 0.05\) и вероятността \(P(\text{Отхвърляне}|H_0=True):= \alpha\) да отхвърлим хипотезата \(H_0\), когато тя е вярна е 0.05 (Грешка от първи род). Областта \(\mathcal{C}\) в която избираме да отхвърлим нулевата хипотеза, когато вектора на наблюденията \((x_1, ..., x_n)\) попадне в нея, наричаме критична област. В нашия случай, критичната област е

\[ \mathcal{C}:= \{(x_1, ..., x_n): \frac{\bar{x}-\mu_{H}}{\sigma/\sqrt{n}} < q^{\mathcal{N}(0,1)}_{\alpha}\} \]

Ако не успеем да отхвърлим хипотезата \(H_0\), а тя се окаже грешна, правиме грешка от втори род. Нейната вероятност \(P(\text{Приемане}|H_0=False)\) бележим с \(\beta\). Сила на теста \(\pi:= 1-\beta\) наричаме вероятността да отхвърлим \(H_0\) при условие, че \(H_0\) наистина е грешна и е вярна \(H_1\).

Вярно Приемаме \(H_0\) Отхвърляме \(H_0\)
\(H_0\) Вярно Грешка от 1ви род с вероятност \(\alpha\)
\(H_1\) Грешка от 2ри род с вероятност \(\beta\) Вярно със сила \(\pi := 1-\beta\)

Когато \(\sigma\) също е неизвестна, то \[ Т := \frac{\bar{X} - 190}{S / \sqrt{n}} \sim \mathcal{T}(n-1) \] и за да видим вероятността да наблюдаваме такава стойност или по-крайна (при услувие, че нулевата хипотеза е вярна), трябва да използваме функцията pt и да заместим \(\sigma\) с \(S\) във формулата за \(Z\)

t <- (mean(heights) - 190) / (sd(heights) / sqrt(100))
pvalue <- pt(t, length(heights) - 1)

# вероятността да сме наблюдавали z или по-крайна стойност от z, при условие, че
# нулевата хипотеза е вярна.
pvalue
## [1] 2.655682e-39

Отново, можем да използваме функцията t.test, която да ни даде цялата информация.

t.test(heights, mu = 190, alternative = "less")
## 
##  One Sample t-test
## 
## data:  heights
## t = -21.451, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is less than 190
## 95 percent confidence interval:
##      -Inf 178.3781
## sample estimates:
## mean of x 
##  177.4031

За разлика от подхода при интервални оценки, при тестване на хипотези, центрираме тестовата статистика \(Z\) около хипотетичен параметър и отхвърляме нулевата хипотеза \(H_0\) когато наблюдаваната (или по-крайна) тест статистика \(Z\) се случва с много малък шанс.

## Picking joint bandwidth of 0.0901

# Тестване на хипотезата, че E(X) < 130
t.test(heights, mu = 130, alternative = "greater")
## 
##  One Sample t-test
## 
## data:  heights
## t = 80.723, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is greater than 130
## 95 percent confidence interval:
##  176.4281      Inf
## sample estimates:
## mean of x 
##  177.4031
## Picking joint bandwidth of 0.0901

# Тестване на хипотезата, че E(X) = 165
t.test(heights, mu = 130, alternative = "two.sided")
## 
##  One Sample t-test
## 
## data:  heights
## t = 80.723, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 130
## 95 percent confidence interval:
##  176.2379 178.5683
## sample estimates:
## mean of x 
##  177.4031

Ако не искаме да правим грешка от първи род, можем да намалим областта на \(\alpha\), но така ще отхвърляме прекалено лесно и ще качим вероятността за грешка от втори род.

Нека симулираме \(T\)-статистиката при вярна нулева хипотеза.

simulate_t <- function(n, th_mean, th_std, hyp_mean) {
  x <- rnorm(n, th_mean, th_std)
  (mean(x) - hyp_mean) / (sd(x) / sqrt(n))
}

# да преброим колко елементи са в критичната област
simulations <- replicate(10000, simulate_t(40, 4, 4, 4))
q <- qt(p = 0.05, df = 39)

mean(simulations <= q)
## [1] 0.0474

Виждаме, че дори при вярна нулева хипотеза, в 1 на 20 случаи ще попадаме в критичната област. Тогава, ако попаднем в тази област при вярна нулева хипотеза, ще отхвърлим нулевата хипотеза като грешна и ще направим грешка от първи род.

Нека симулираме данни от две различни разпределение ; \(X_1 \sim \cal{N(\mu_1, \sigma^2)}\) и \(X_2 \sim \cal{N(\mu_2, \sigma^2)}\) с цел да видим как се променят грешките от втори род \(\beta\) при увеличаване на броя на елементите \(n\) с които изграждаме тестовата статистика.

simulate_t <- function (n, m, std = 10) {
  observations <- rnorm(n, m, std)
  mean(observations) / (sd(observations) / sqrt(n))
}

simulate_ht <- function(n, n_sim, mean1, mean2, std = 10) {
  set.seed(1)

  bind_rows(
    data.frame(
      x = replicate(n_sim, simulate_t(n, mean1, std)),
      distribution = paste0("средно=", mean1),
      n = as.factor(n)
    ),
    data.frame(
      x = replicate(n_sim, simulate_t(n, mean2, std)),
      distribution = paste0("средно=", mean2),
      n = as.factor(n)
    )
  )
}

c(5, 10, 20, 30, 70, 140) %>%
  imap(~simulate_ht(.x, 10000, 5, 8, 6)) %>%
  bind_rows() %>%
  ggplot(mapping = aes(x = x, fill = distribution)) +
  geom_density(alpha = 0.55) +
  xlim(c(-5, 20)) +
  facet_wrap(~n) +
  xlab("Т-статистика") +
  ylab("Вероятностни плътности") +
  ggtitle("Промяна на силата при увеличаване на броя на елементи") +
  theme_bw() +
  theme(legend.position = "bottom")
## Warning: Removed 48 rows containing non-finite values (stat_density).

Виждаме, че при растящ брой елементи \(n\), областта в която бихме могли да направим грешка от втори род намалява, понеже дисперсията на T-статистиката под нулевата хипотеза намалява и разпределението се центрира около хипотетичния си параметър. По този начин, \(\beta\) намалява, докато силата на теста \(\pi:= 1-\beta\) се увеличава.

Забелязваме, че силата на теста \(\pi\) расте при увеличаване на броя на наблюденията, но при по-голямо стандартно отклонение, скоростта на растеж на \(\pi\) намалява.

При качване на нивото на значимост \(\alpha\), качваме силата на теста (\(\pi\)) и намаляваме вероятността за грешка от втори род \(\beta\), но по този начин увеличаваме вероятността за грешка от първи род \(\alpha\).

Тестване на хипотези върху две извадки

Ако имаме еднакъв брой наблюдения от две случайни величини \(X\) и \(Y\), които са нормално разпределени с еднакви дисперсии - съответно, \(X \sim \mathcal{N}(\mu_1, \sigma^2)\) и \(Y \sim \mathcal{N}(\mu_2, \sigma^2)\), можем да зададем нулева хипотеза \(H_0 := \mu_1 = \mu_2\). Ако нулевата хипотеза е вярна, то \(\bar{X} - \bar{Y} \sim \mathcal{N}(0, \frac{2\sigma^2}{n})\). Тогава

\[ Т := \frac{\bar{X}-\bar{Y}}{\sqrt{\frac{S_1^2}{n} + \frac{S_2^2}{n}}} \sim \mathcal{T}(2n - 2) \]

Тогава, можем да използваме досегашния подход при тестване на хипотези за да проверим дали наистина \(\mu_1=\mu_2\).

Aко имаме \(pvalue < \alpha\), където \(\alpha\) е нивото на значимост, а t e наблюдаваната стойност за T-статистиката, можем да отхвърлим нулевата хипотеза, че \(\mu_1=\mu_2\).

Пример:
# брой елементи в двете случайни извадки
n <- 40

# Симулации на случайни извадки от с.в. X ~ N(3, 9) и Y ~ N(4, 9)
x <- rnorm(n, 3, 3)
y <- rnorm(n, 4, 3)

# Извадкови вариации
v1 <- var(x)
v2 <- var(y)

# Степени на свобода за T-разпределението
degrees_of_freedom <- 2 * n - 2

# Наблюдаваната стойност от T-статистиката
t <- (mean(x) - mean(y)) / sqrt(v1 / 40 + v2 / 40)

# Изчисление за pvalue при двустранен тест
pval <- 2 * min(
  pt(q = t, df = degrees_of_freedom, lower.tail = TRUE),
  pt(q = t, df = degrees_of_freedom, lower.tail = FALSE)
)

# резултат
c("t" = t, "pvalue" = pval, "df" = degrees_of_freedom)
##          t     pvalue         df 
## -1.9004828  0.0610646 78.0000000

Или ако просто използваме t.test

t.test(x, y, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  x and y
## t = -1.9005, df = 78, p-value = 0.06106
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.67067950  0.06201836
## sample estimates:
## mean of x mean of y 
##  2.835611  4.139941

Когато работим с предположението, че дисперсиите са еднакви, но имаме различен брой наблюдения от двете случайни извадки, за да оценим по-точно общата дисперсия посредством извадъчните дисперсии, използваме формулата за комбинирана извадъчна дисперсия и съответно - фoрмулата за комбинирано стандартно отклонение.

\[ S_p^2 = \frac{(n - 1)S_1^2 + (m - 1)S_2^2}{n + m - 2} \]

където \(n\) и \(m\) са броят на елементите в двете извадки, а \(S_i\) са извадъчните дисперсии на двете извадки. В този случай T-статистиката е

\[ Т := \frac{\bar{X}-\bar{Y}}{S_p\sqrt{\frac{1}{n} + \frac{1}{m}}} \sim \mathcal{T}(n + m - 2) \]

Пример:
# Формула за смесено извадъчно стандартно отклонение
pooled_sd <- function (x, y) {
  n1 <- length(x)
  n2 <- length(y)
  sqrt(
    ((n1 - 1) * var(x) + (n2 - 1) * var(y))
    /
    (n1 + n2 - 2)
  )
}

# броя на елементите в двете случайни извадки които ще симулираме
n <- 30
m <- 40

# Симулации с различен брой над с.в.
x <- rnorm(n, 3, 2)
y <- rnorm(m, 4, 2)

# Степени на свобода
degrees_of_freedom <- n + m - 2

# Наблюдавана стойност t на Т-статистиката
t <- (mean(x) - mean(y)) / (pooled_sd(x, y) * sqrt((1 / n) + (1 / m)))

# Изчисление за pvalue при двустранен тест
pvalue <- 2 * min(
  pt(t, degrees_of_freedom, lower.tail = TRUE),
  pt(t, degrees_of_freedom, lower.tail = FALSE)
)

# резултат
c("t" = t, "pvalue" = pvalue, "df" = degrees_of_freedom)
##           t      pvalue          df 
## -1.86249033  0.06685454 68.00000000

Или с t.test

t.test(x, y, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  x and y
## t = -1.8625, df = 68, p-value = 0.06685
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.80818636  0.06232572
## sample estimates:
## mean of x mean of y 
##  2.937953  3.810884

Когато работим с предположението, че дисперсиите на двете случайни величини \(X\) и \(Y\) не са еднакви, тоест \(\sigma_X^2 := D(X) \neq D(Y) =: \sigma_Y^2\), за степените на свобода използваме формулата за комбинирани степени на свобода.

\[ df_{c} := \frac{(\frac{S_1^2}{n} + \frac{S_2^2}{m})^2}{\frac{(S_1^2/n)^2}{n-1} + \frac{(S_2^2/m)^2}{m-1}} \]

В този случай Т-статистиката ще е

\[ Т := \frac{\bar{X}-\bar{Y}}{\sqrt{\frac{S_1^2}{n} + \frac{S_2^2}{m}}} \sim \mathcal{T}(df_{c}) \]

Пример:
# брой наблюдения над симулираните случайни извадки
n <- 30
m <- 40

# Симулации с различен брой наблюдения над нормално разпределени с.в.
x <- rnorm(n, 3, 2)
y <- rnorm(m, 4, 3)

# Извадъчни дисперсии
v1 <- var(x)
v2 <- var(y)

# Смесени степени на свобода
degrees_of_freedom <- (
  (v1 / n + v2 / m)^2
  /
  ((v1 / n)^2 / (n - 1) + (v2 / m)^2 / (m - 1)))

t <- (mean(x) - mean(y)) / sqrt(v1 / n + v2 / m)

# Изчисление за pvalue при двустранен тест
pvalue <- 2 * min(
  pt(t, degrees_of_freedom, lower.tail = TRUE),
  pt(t, degrees_of_freedom, lower.tail = FALSE)
)

# резултат
c("t" = t, "pvalue" = pvalue, "df" = degrees_of_freedom)
##         t    pvalue        df 
## -1.151663  0.253680 64.971366

Или за по-кратко с t.test

t.test(x, y, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = -1.1517, df = 64.971, p-value = 0.2537
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.9534532  0.5245233
## sample estimates:
## mean of x mean of y 
##  3.077375  3.791840

Сдвоени и несдвоени данни (Paired vs Unpaired)

Сдвоени данни (наблюдения) от две извадки наричаме, такива при при които имаме зависимост между извадките и няблюденията идват от едни и същи индивиди, но в различен етап от време. Тоест \(x_i\) и \(y_i\) са наблюдения върху един и обект. Пример за това би бил ефекта върху нивото на холестерола преди и след приемането на лекарство от едни и същи хора. При несдвоени данни, имаме наблюдения от независими извадки \(X\) и \(Y\). Пример за това би бил, кръвното на група която не приема определено лекарство, и кръвното на група която е приела лекарството. При t.test в този случай, при нулева хипотеза \(H_0 := \mu_0=0\), тест статистиката \(Т\) би била \[ T := \frac{(\bar{X} - \bar{Y})-(\mu_1-\mu_2)}{st.dev(\bar{X}-\bar{Y})} = \frac{\bar{X}_D-\mu_0}{S_D/\sqrt{n}} \sim \mathcal{T}(n - 1) \] където \(\bar{X}_D\) и \(S_D\) са извадъчното средно и стандартно отклонение върху разликите между всички наблюдения.

Пример:
# Брой елементи от които ще симулираме извадка
n <- 30

# Симулираме наблюдения от с.в.
x <- rnorm(n, 80, 4)
y <- rnorm(n, 84, 4)

# Тест статистика при сдвоен тест
t <- mean(x - y) / (sd(x - y) / sqrt(30))

# Изчисление на pvalue
pvalue <- 2 * min(pt(t, n - 1), pt(t, n - 1, lower.tail = FALSE))

# резултат
c("t" = t, "pvalue" = pvalue, "df" = n - 1)
##           t      pvalue          df 
## -2.38593558  0.02378466 29.00000000

Отново, можем да ползваме t.test

t.test(x, y, paired = TRUE)
## 
##  Paired t-test
## 
## data:  x and y
## t = -2.3859, df = 29, p-value = 0.02378
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.5515669 -0.3499634
## sample estimates:
## mean of the differences 
##               -2.450765

Кога ползваме t.test? Тест на Уилкоксън

Т-тестът работи с предположението, че имаме извадки от нормално разпределени случайни величини (популации) и в този случай, трябва да проверим дали нашите данни отговарят на това условие. Поради Централната Гранична Теорема, стига да имаме достатъчно наблюдения (като правило се счита над 30), можем да ползваме Т-тест и без данните ни да следват нормално разпределение, защото

\[ \frac{\bar{X} - \mu}{\sigma/\sqrt{n}} \to_d \mathcal{N}(0, 1) \]

Когато тези условия не са спазени, можем да ползваме Непараметричен тест на Уилкоксън.

Тест на Уилкоксън

Тестът на Уилкоксън е непараметричен тест, който ни позволява да да сравним две извадки, без да правим преположение за вероятностното разпределение което стои зад тях. Предположенията върху които той стъпва са следните

  • Наблюденията са независими едно от друго (Случайна Извадка)
  • Наблюденията са числени или поне ординални (можем да ги сравним).

При теста на Уилкоксън, работим с нулевата хипотеза \(H_0 := P(x_i > y_i) = \frac{1}{2}\), срещу алтернативата \(H_1 := P(x_i > y_i) \neq \frac{1}{2}\) и следната \(U\)-статистика.

Нека имаме \(X_1, ...,X_n\) е независима случайна извадка от \(X\) и \(Y_1, ...,Y_n\) е независима случайна извадка от \(Y\) и

\[ U := \sum_{i=1}^n\sum_{j=1}^mS(X_i, Y_j) \]

където

\[ S(X, Y) := \begin{cases} 1,& \text{if } Y < X\\ 1/2, & \text{if } Y = X \\ 0, & \text{if } Y > X \\ \end{cases} \] Разпределението на \(U\) е добре известно и го има във всеки справочник. В R функцията на разпределението му е pwilcox. Виждаме, че \(U\) брои броя на инверсиите между Y и Х. Тоест, ако имаме изместеност (Y и X не са с един център), наблюдаваната стойност за \(U\)-статистиката ще бъде в някаква крайна област на разпределението.

Пример:
s <- function(x, y) (y < x) + .5 * (y == x)

x <- c(1, 5, 7,  2, 5,  9, 3, 2, 1, 12)
y <- c(8, 3, 13, 11, 6, 12, 15, 10, 6, 8)

U <- sum(unlist(map(x, ~sum(s(.x, y)))))

c("U" = U, "pvalue" = 2 * pwilcox(U, n = length(x), m = length(y)))
##           U      pvalue 
## 18.00000000  0.01468964

Отново, има функция която изпълнява целия тест

# Дава леко по-различни резултати поради вътрешни оптимизации
wilcox.test(x, y)
## Warning in wilcox.test.default(x, y): cannot compute exact p-value with ties
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  x and y
## W = 18, p-value = 0.01696
## alternative hypothesis: true location shift is not equal to 0

Виждаме, че разпределението прилича на Нормално разпределение при увеличение на броя на наблюдения \(n\). При достатъчен брой наблюдения, разпределението на тест статистиката може да се апроксимира с нормално.

Тестът на Уилкоксън може да бъде изпълнен при всички от горните случаи

  • само върху една извадка - wilcox.test(x, mu = 10, alternative = “two.sided”)
  • върху две извадки - wilcox.test(x, y, alternative = “two.sided”)
  • сдвоени данни - wilcox.test(x, y, alternative = “two.sided”, paired = TRUE)

\(\chi^2\)-тест на Пиърсън

Интуиция за две групи

Ако имаме случайна извадка \((Y_1, ..., Y_n)\) от биномно разпределена сл.в. \(Y \sim Bin(n, p)\), можем да обозначим броя наблюдавани успехи (елементи от едната група) с \(О_S\), а броя наблюдавани неуспехи с \(О_F\) (елементи от другата група), и съответно техните очаквани стойности с \(E_S\) и \(E_F\).

Тогава, от Централната гранична теорема

\[ X := \frac{Y-np}{\sqrt{np(1-p)}} \to_d \mathcal{N}(0, 1) \]

Следователно \(X^2 \to_d \chi^2(1)\) и ако имаме достатъчно наблюдения \(n\), \(X^2\) ще има приблизително \(\chi^2\) разпределение с една степен на свобода.

Oт това, че \(\frac{1}{p(1-p)}=\frac{1}{p} + \frac{1}{n-p}\), и това, че \((Y-np)^2=[(n-Y)-n(1-p)]^2\), можем да представим \(X^2\) по следния начин \[ X^2 := \frac{(Y-np)^2}{np(1-p)} = \frac{(Y-np)^2}{np}+\frac{(Y-np)^2}{n(1-p)}=\\ = \frac{(Y-np)^2}{np}+\frac{[(n-Y)-n(1-p)]^2}{n(1-p)} =\\ = \frac{(O_S-E_S)^2}{E_S}+\frac{(O_F-E_F)^2}{E_F} \] По този начин, при \(n\) наблюдавани стойности, можем да тестваме нулева хипотеза, че очакваният брой успехи (елементи от едната група) \(E_S = np_S\) е определено число, а очаквания брой неуспехи ( елементи от другата група) \(E_F =np_F\)е друго число. Също, можем да тестваме и хипотези за самите вероятности на двете групи \(p_S\) и \(p_F\). Ако нулевата хипотеза не е вярна, наблюдаваната \(X^2\) статистика ще е прекалено надясно от \(0\).

Както при биномно разпределена случайна величина, така и при такава с много групи (мултиномно разпределена), съществува асимптотична апроксимация с нормално разпределение. В такъв случай, можем да направим сходни разсъждения за \(k\) на брой различни групи и \[ X^2 := \sum_{i=1}^k\frac{(O_i-E_i)^2}{E_i} \]

ще има, отново асимптотично приблизително \(\chi^2\) разпределение, но този път с \(k-1\) степени на свобода, където \(О_i\) и \(E_i\) са наблюдавания и очаквания брой елементи под \(H_0\) от \(i\)-тата група.

Група 1 Група 2 … Група k
Наблюдавани \(O_1 := n\hat{p_1}\) \(O_2 := n\hat{p_2}\) … \(O_k := n\hat{p_k}\)
Очаквани под \(H_0\) \(E_1 := np_1\) \(E_2 := np_2\) … \(E_k := np_k\)

Нека симулираме разпределенията на \(X^2\) при вярна нулева хипотеза и различен брой групи, с цел да видим, че наистина имат приблизително \(\chi^2\) разпределения.

simulate_chsq <- function(n = 60, th_probs, hyp_probs) {
  k = length(hyp_probs)
  # Правим извадка от к на брой групи с теоретични вероятности за всяка - th_probs
  x <- sample(1:k, prob = th_probs, size = n, replace = TRUE)

  # изграждаме брой елементи във всяка група
  observed <- unlist(map(1:k, ~sum(x == .x)))

  # очакваните бройки елементи ако нулевата хипотеза е вярна
  expected <- n * hyp_probs

  # връщаме X^2 статистиката
  sum((observed - expected)^2 / expected)
}

# общ брой елементи от всички групи
n <- 120
# Брой симулации на X^2-статистиката
s <- 2000

# за степени на свобода от 4 до 10 генерираме X^2 статистики
5:11 %>%
  imap(
    ~data.frame(
      # Задаваме истинските вероятности за нулева хипотеза, така че да видим разпределенията
      x = replicate(s, simulate_chsq(n, th_probs = rep(1 / .x, .x), hyp_probs = rep(1 / .x, .x))),
      df = as.factor(.x)
    )
  ) %>%
  # сливаме всички data.frames в един data.frame с цел да визуализираме по-лесно
  bind_rows() %>%
  ggplot(mapping = aes(x = x, fill = df)) +
  geom_density(alpha = 0.25) +
  ggtitle("X^2 статистиката с Хи-квадрат разпределение") +
  xlab("X^2") +
  ylab("Емпирична Плътност")

Нека направим тест с четири групи симулирани с теоретични вероятности 1/4 за всяка и \(n=60\). Тогава, ако тестваме нулевата хипотеза \(H_0 := \{p_1=1/2, p_2 = p_3=p_4=1/6\}\), нашата \(X^2\) статистика би трябвало да е прекалено надясно от нулата.

xsquared <- simulate_chsq(n = 60, th_probs = rep(1 / 4, 4), hyp_probs = c(1/2, 1/6, 1/6, 1/6))

pvalue <- pchisq(xsquared, df = 3, lower.tail = FALSE)

c("X^2" = xsquared, "pvalue" = pvalue, "df" = 3)
##          X^2       pvalue           df 
## 15.800000000  0.001246226  3.000000000
Група 1 Група 2 Група 3 Група 4
Истински Очаквания 15 15 15 15
\(O_i := n\hat{p_i}\) := Наблюдавани 14 13 17 16
\(E_i := np_i\) := Очаквани под \(H_0\) 30 10 10 10
\(\frac{(O_i-E_i)^2}{E_i}\) \(\frac{(14-30)^2}{30}=8.53\) \(\frac{(13-10)^2}{10}=0.90\) \(\frac{(17-10)^2}{10}=4.90\) \(\frac{(16-10)^2}{10}=3.60\)

Виждаме, че вероятността да наблюдаваме стойност по-голяма или равна от \(17.93\) на \(X^2\)-статистиката е много малка (\(pvalue := P(X^2\ge17.93) < 0.005\)). В такъв случай, можем да отхвълим нулевата хипотеза, че очакваните вероятности са \(\{p_1=1/2, p_2 = p_3=p_4=1/6\}\).

Критичните области и техните pvalue-та биха изглеждали така

## Picking joint bandwidth of 1.71

В R за улеснение, използваме функцията chisq.test(x, p) където x ще е наблюдаван брой елементи в различните групи, а p е вероятностите под нулевата хипотеза, за различните групи.

Приложения на \(\chi^2\)-тест

Проверка за произволно дискретно разпределение

Ако имаме дискретни данни, бихме могли да тестваме хипотезата, че наблюдаваната извадка идва от дискретна случайна величина с избрано от нас разпределение.

Пример

Нека имаме случайна извадка \(X_1, ..., X_{100} \sim^d X \sim Poi(2)\). Искаме да тестваме хипотезата \(H_0 := \lambda = 2.5\) срещу алтернативата \(H_1 := \lambda \neq 2.5\)

# брой наблюдения които имаме
n <- 100

set.seed(2)
# Стойностите които примерно сме наблюдавали
observed_props <- rpois(n, lambda = 2) %>% 
  table() %>%
  prop.table() %>%
  unname() %>%
  as.vector()

# Хипотетични вероятности под нулева хипотеза
hypothet_props <- dpois(0:6, lambda = 2.5)

# Калибрация с цел по-лесно тестване. chisq.test изисква вероятностите да се сумират до 1
hypothet_props <- hypothet_props + (1 - sum(hypothet_props)) / 7

# Визуализация на пропорциите
bind_rows(
  data.frame(x = 0:6, y = hypothet_props, type = "Хипотетични"),
  data.frame(x = 0:6, y = observed_props, type = "Наблюдавани")
)  %>%
  ggplot(mapping = aes(x = x, y = y, fill = type)) +
  geom_bar(position = "dodge", stat = "identity") +
  ggtitle("Наблюдавани и хипотетични вероятности за разпределението") +
  xlab("X") +
  ylab("Плътност") +
  theme_light()

Така, ако нулевата хипотеза е вярна,

\[ X^2 := \sum_{i=0}^6\frac{(n\hat{p_i}-np_i)^2}{np_i} \sim^d_{aprox} \chi^2_{6} \]

observed <- n * observed_props
expected <- n * hypothet_props

xsquared <- sum((observed - expected)^2 / expected)
pvalue <- pchisq(xsquared, lower.tail = FALSE, df = 6)

c("X^2" = xsquared, "df" = 6, "pvalue" = round(pvalue, 6))
##       X^2        df    pvalue 
## 16.964442  6.000000  0.009415

Така виждаме, че е много малко вероятно да наблюдаваме съответната тестова статистика, ако истинското разпоределение беше \(Poi(2.5)\). В такъв случай, бихме могли да отхвърлим нулевата хипотеза при избрано от нас ниво на съгласие \(\alpha\).

В R за по-лесно това би станало така

chisq.test(observed, p = hypothet_props)
## Warning in chisq.test(observed, p = hypothet_props): Chi-squared approximation
## may be incorrect
## 
##  Chi-squared test for given probabilities
## 
## data:  observed
## X-squared = 16.964, df = 6, p-value = 0.009415

Проверка за произволно непрекъснато разпределение

Ако имаме непрекъснати данни и искаме да тестваме хипотезата, че те идват от случайна величина с дадено непрекъснато разпределение, бихме могли да дискретизираме данните в под интервали и да проверим дали честотата с която се падат е същата като при същата дискретизация на плътността на тестваното от нас разпределение.

Пример

Нека имаме \(2000\) на брой наблюдения \(X_1, ..., X_{2000}\) от случайна величина \(X \sim \Gamma(12, 4)\) и нека тестваме нулева хипотеза с истинските параметри \(\alpha = 12, \beta = 4\). В този случай, не би трябвало да можем да отхвърлим нулевата хипотеза.

# Брой наблюдения
n <- 2000

# Параметри под нулевата хипотеза за разпределението
hyp_shape <- 12; hyp_rate <- 4;

# seed с цел възпроизвеждане на резултатите
set.seed(16)

# Симулирани данни от Гама разпределена случайна величина
x <- rgamma(n, 12, 4)

# Дискретизация
divisions <- seq(0.5, 6.5, by = 0.5)

# Хипотетични честоти ако нулевата хипотеза е вярна
hyp_props <- pgamma(divisions, hyp_shape, hyp_rate) - pgamma(divisions - 0.5, hyp_shape, hyp_rate)

# таблица с брой наблюдения в интервалите
observed_table <- x %>%
  cut(c(0, divisions)) %>%
  table()

# Наблюдавани честоти
observed <- as.vector(unname(observed_table))

# Визуализация на честотите (подобно е на хистограма)
bind_rows(
  data.frame(x = names(observed_table), y = observed, type = "Наблюдавани"),
  data.frame(x = names(observed_table), y = n * hyp_props, type = "Хипотетични"),
) %>%
  ggplot(mapping = aes(x = x, y = y, fill = type)) +
  geom_bar(position = "stack", stat = "identity") +
  coord_flip() +
  xlab("Интервали на Дискретизация") +
  ylab("Честота на срещане в интервала") +
  ggtitle("Честота на Дискретизирани данни")

# добавяме остатъка за да се сумират до 1
hyp_props <- hyp_props + (1 - sum(hyp_props)) / length(hyp_props)

# изпълнение на теста
chisq.test(observed, p = hyp_props)
## Warning in chisq.test(observed, p = hyp_props): Chi-squared approximation may be
## incorrect
## 
##  Chi-squared test for given probabilities
## 
## data:  observed
## X-squared = 7.6662, df = 12, p-value = 0.8106

Не успяваме да отхвърлим нулевата хипотеза, че данните са наблюдения от случайна величина \(X \sim \Gamma (12, 4)\).

Проверка за независимост на две случайни величини

Ако имаме две случайни извадки \(X_1, ..., X_n \sim^d X\) и \(Y_1, ..., Y_n \sim^d Y\), бихме могли отново да изградим таблици с пропорции \(\hat{p}_{ij} :=\hat{P}(X = i, Y = j)\) от наблюдаваните стойности, и съответно маргиналните пропорции (емпирични вероятности)
\(\hat{p}_{i, \forall}:= \hat{P}(X=i)=\sum_{j = 1}^c\hat{p}_{ij}\) и \(\hat{p}_{\forall, j}:= \hat{P}(Y=j)=\sum_{i = 1}^r\hat{p}_{ij}\), където \(r\) и \(c\) са бройките елементи за \(X\) и \(Y\).

В такъв случай, ако случайните величини \(X\) и \(Y\) са независими, то \(P(X = x_i, Y = y_j) = P(X = x_i)\cdot P(Y=y_j) := p_{i, \forall}p_{\forall,j}\) за всяко \(i, j\) и би трябвало да наблюдаваме емпирични вероятности които приблизително притежават това свойство.

Тогава

\[ X^2 := \sum_{i=1}^r\sum_{j=1}^c \frac{(O_{i,j}-E_{i,j})^2}{E_{i,j}} \sim^d_{approx} \chi^2_{(r-1)(c-1)} \]

Пример

Нека симулираме две извадки с големина \(n=600\) от независими случайни величини \(X \sim Bin(8, 0.5), Y \sim Poi(2.3)\) и изпълним \(\chi^2\)-тест за независимост, като евентуално след това, бихме могли да добавим зависимост между част от наблюденията с цел да видим как се променят резултатите от \(\chi^2\) теста.

# генериране на две независими случайни величини
x <- rbinom(600, size = 8, prob = 1 / 2)
y <- rpois(600, 2.3)

# Вкарване на зависимост
y[400:600] <- y[400:600] + x[400:600]

# таблица брояща наблюдавани комбинации от вида (X, Y)
observed <- table(x, y)

# брой на всички елементи, N = n * m
N <- sum(observed)

# Маргинални емпирични вероятности
x_marginal <- prop.table(table(x))
y_marginal <- prop.table(table(y))

# Очаквани бройки при предположение за независимост
expected <- N * (x_marginal %*% t(y_marginal))

# X^2 статистика
xsquare <- sum((observed - expected)^2 / expected)

# Степени на свобода df = (r - 1) * (s - 1)
degrees_of_freedom <- prod(dim(observed) - 1)

# pvalue за съответната наблюдавана Х^2 статистика
pvalue <- pchisq(xsquare, df = degrees_of_freedom, lower.tail = FALSE)

# връщане на резултат
c("X^2" = xsquare, "df" = degrees_of_freedom, "pvalue" = pvalue)
##          X^2           df       pvalue 
## 1.461809e+02 1.040000e+02 4.054662e-03

Или за кратко с функцията chisq.test

chisq.test(observed)
## Warning in chisq.test(observed): Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  observed
## X-squared = 146.18, df = 104, p-value = 0.004055

Функцията chisq.test автоматично разбира, че става въпрос за тест за независимост, когато ѝ се подаде таблица.

Проста Линейна Регресия

Често се налага да обясним една непрекъсната променлива чрез друга, примерно

Нека разгледаме данни от последния пример

Виждаме, че има добра линейна връзка между печалбите и приходите вложени в ТВ реклами. Данните като цяло се движат около права линия около която има лек случаен шум. Ако моделираме печалбата като случайна величини \(Y\) и \(X\), бихме могли да кажем, че

\[ Y = \beta_0 + \beta_1X + \epsilon, \] където \(\epsilon\) е нашия случаен елемент, който за леснота можем да считаме като независими наблюдения oт нормално разпределена случайна величина, тоест

\[ \epsilon_1, ..., \epsilon_n \sim^d \epsilon \sim \mathcal{N}(0, \sigma^2) \]

За леснота, можем да считаме, че \(Х\) са наблюдавани константи и техните вероятностни свойства не ни интересуват.

По този начин,

\[ E[Y] = E(\beta_0 + \beta_1x + \epsilon) = \beta_0 + \beta_1x\\ Var[Y] = Var(\beta_0 + \beta_1x + \epsilon) = \sigma^2\\ Y \sim \mathcal{N}(\beta_0+\beta_1x, \sigma^2) \]

и математическото очакването за \(Y\) ще е права линия около която наблюденията варират.

Един начин по който бихме могли да оценим математическото очакване на \(Y\) е метода на най-малките квадрати.

Понеже знаем, че върху математическото очакване може да се гледа като на център на тежестта, естествен начин за неговата оценка е да изберем минималното разстояние между наблюденията \(\overrightarrow{y} = (y_1, ..., y_n)'\) и линията \(\beta_0 + \beta_1x\), тоест за \(\beta_0\) и \(\beta_1\) да изберем такива коефициенти, че

\[ (\hat{\beta_0}, \hat{\beta_1}) := \underset{\beta_0, \beta_1}{argmin}\text{ }\rho(\overrightarrow{y}, \beta_0 + \beta_1 \overrightarrow{x})^2=\\ =\underset{\beta_0, \beta_1}{argmin}\text{ }\left\lVert \overrightarrow{y}-\beta_0 - \beta_1 \overrightarrow{x} \right\rVert^2_2=\\ =\underset{\beta_0, \beta_1}{argmin}\text{ }\sum_{i=1}^n(y_i-\beta_0-\beta_1 x_i)^2 \]

Тоест, избираме коефициенти \(\beta_0\) и \(\beta_1\), такива че нашата линия \(f(x) := \beta_0 + \beta_1x\) максимално да приближава наблюденията \((y_1, ..., y_n)\) които имаме. Функцията която минимизираме по-горе, наричаме сума на квадратите на грешките и най-често обозначаваме с RSS,

която като функция на \(\beta_0\) и \(\beta_1\) би изглеждала така

Един начин по който можем да намерим минимума на горната функция е като намерим първите частни производни и ги положим на 0, т.е.

\[ \frac{\partial{(RSS)}}{\beta_0}=0\\ \frac{\partial{(RSS)}}{\beta_1}=0 \]

и да намерим инфлексната точка на функцията. Решавайки за \(\beta_0\) и \(\beta_1\) получаваме

\[ \hat{\beta_0} = \bar{y}-\beta_1\bar{x}\\ \]

и съответно

\[ \hat{\beta_1} = \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i-\bar{x})^2}\\ \]

Забелязваме, че оценката на \(\beta_1\) може да се представи и така

\[ \hat{\beta_1} = \langle\frac{(\overrightarrow{x}-\bar{x})}{||\overrightarrow{x}-\bar{x}||_2},\frac{(\overrightarrow{y}-\bar{y})}{||\overrightarrow{x}-\bar{x}||_2} \rangle \]

Тоест, оценката за \(\beta_1\) ни дава колко наблюдаваните вектори \(\overrightarrow{x}\) и \(\overrightarrow{y}\) (центрирани и нормирани спрямо \(\overrightarrow{x}\)) са в една посока, или колко \(\overrightarrow{y}\) се променя спрямо \(\overrightarrow{x}\) (скаларно произведение).

По този начин, получената оценка за \(EY\) e \[\hat{Y} :=: \hat{f}(x) := \hat{\beta_0}+\hat{\beta_1}x\]

## `geom_smooth()` using formula 'y ~ x'

Свойства на \(\hat{\beta}_0\) и \(\hat{\beta}_1\)

Нека разгледаме свойствата на оценките \(\hat{\beta_0}\) и \(\hat{\beta_1}\).

\[ Е[\hat{\beta}_1] = E[\frac{\sum_{i=1}^n(x_i-\bar{x})(Y_i-\bar{Y})}{\sum_{i=1}^n(x_i-\bar{x})^2}] =\\ =\frac{\sum_{i=1}^n(x_i-\bar{x})E[Y_i]}{\sum_{i=1}^n(x_i-\bar{x})^2}=\\ =\frac{\sum_{i=1}^n(x_i-\bar{x})(\beta_0 + \beta_1x_i)}{\sum_{i=1}^n(x_i-\bar{x})^2}=\\ =\frac{\sum_{i=1}^n(x_i-\bar{x})\beta_0}{\sum_{i=1}^n(x_i-\bar{x})^2}+\frac{\beta_1\sum_{i=1}^n(x_i-\bar{x})(x_i-\bar{x}+\bar{x})}{\sum_{i=1}^n(x_i-\bar{x})^2}=\\ 0+\frac{\beta_1\sum_{i=1}^n(x_i-\bar{x})^2}{\sum_{i=1}^n(x_i-\bar{x})^2}+0=\\ =\beta_1 \]

и така, оценката \(\hat{\beta_1}\) за \(\beta_1\) е неизместена оценка.

\[ Е[\hat{\beta_0}]=E[\bar{y}-\hat{\beta_1}\bar{x}]=E[\beta_0+\beta_1\bar{x}+\bar{\epsilon}-\hat{\beta_1}\bar{x}]=\beta_0 \]

и следователно \(\hat{\beta}_0\) също е неизместена оценка за \(\beta_0\).

\[ Var[\hat{\beta}_1] = Var[\frac{\sum_{i=1}^n(x_i-\bar{x})(Y_i-\bar{Y})}{\sum_{i=1}^n(x_i-\bar{x})^2}]=\\ =\frac{\sum_{i=1}^n(x_i-\bar{x})^2\sigma^2}{[\sum_{i=1}^n(x_i-\bar{x})^2]^2}=\frac{\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2} \]

За дисперсията на \(\hat{\beta}_0\) получаваме

\[ Var[\hat{\beta}_0]=Var[\bar{y}-\hat{\beta_1}\bar{x}]=\\ =\frac{\sigma^2}{n}+\bar{x}^2Var[\hat{\beta_1}]+Cov[\bar{y}, \hat{\beta_1}]=\\ =\frac{\sigma^2}{n}+\frac{\bar{x}\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2} \]

понеже \(Cov[\bar{y}, \hat{\beta_1}] = 0\) (Проверете сами!).

За ковариацията между \(\hat{\beta}_0\) и \(\hat{\beta}_1\) получаваме

\[ Cov(\hat{\beta}_0, \hat{\beta}_1) = E[(\bar{y}-\hat{\beta}_1\bar{x})\hat{\beta}_1] - \beta_0\beta_1= E[\bar{y}]E[\hat{\beta}_1]-\bar{x}E[\hat{\beta}_1^2]=\\ =(\beta_0+\beta_1\bar{x})\beta_1-\bar{x}Var[\hat{\beta}_1]-\bar{x}\beta_1^2-\beta_0\beta_1=\\ =-\bar{x}Var[\hat{\beta}_1]=-\frac{\bar{x}\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2} \]

Оценка за \(\sigma^2\)

Нека разгледаме минимизираната стойност на \(RSS\) като случайна величина

\[ \hat{RSS} := ||\overrightarrow{Y}- \overrightarrow{\hat{Y}}||_2^2 = \sum_{i=1}^n{(Y_i-\hat{Y}_i)^2} \]

Тогава

\[ \frac{\hat{RSS}}{\sigma^2} \sim \chi^2_{n-2} \] и ако вземем математическото очакване на \(\hat{RSS}\)

\[ E[\frac{\hat{RSS}}{\sigma^2}] = n-2 \implies E[\frac{\hat{RSS}}{n-2}] = \sigma^2 \]

Така получаваме неизместена оценка \(\hat{\sigma}^2 := \frac{\hat{RSS}}{n-2}\) с дисперсия \(2\sigma^4/(n-2)\).

Интервали на достоверност в Проста Линейна Регресия

Тъй като стойността на \(\hat{y}_i\) във всяка точка \(x_i\) e оценка за математическото очакване в тази точка, бихме могли да изградим интервал на достоверност

## `geom_smooth()` using formula 'y ~ x'